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Testing  finite  differencing  schemes  for 
the  shallow  water  equations 

G.  Peggion 


Abstract:  The  explicit,  the  semi-implicit,  and  the  fractional  step  sche¬ 
mes  are  tested  and  compared  in  the  solution  of  the  shallow  water  equations 
The  explicit  finite-difference  formulation  is  the  most  accurate,  but  is  restric¬ 
ted  by  a  stability  condition  which  is  not  suitable  for  long-term  numerical 
simulations.  The  standard  semi-implicit  scheme  requires  the  solution  of  an 
elliptic  equation  which  is  also  time-consuming.  The  fractional  step  method 
results  in  the  least  accurate,  but  computationally  the  most  efficient  solu 
tion. 


Keywords:  Courant-Friedrichs-Levy  stability  condition  •  explicit 

scheme  ©  fractional  step  method  o  semi-implicit  scheme  •  shallow 
water  equation 
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1.  Introduction 


I  he  study  of  oo-aii  and  atmosphere  dynamics  classically  develops  from  setting 
mathematical  frameworks.  These  may  he  collections  <>f  numhers  labelled  as  obser¬ 
vations  or  systems  of  equations,  considered  as  ‘truly  representative’  of  the  dynamics 
Within  the  limits  of  the  assumptions  and  approximations  of  the  analytical  model, 
the  numerical  technique  used  for  solving  the  problem  may  introduce  two  additio 
nal  distortions  into  the  representation  of  the  solutions:  the  error  inherent  from  the 
truncated  arithmetic  and  the  error  created  by  approximating  continuous  differential 
equations  with  discrete  algebraic  expressions  (Cirotjahn  and  O'Brien,  107**)  Howe 
ver.  the  accuracy  of  the  solutions  is  not  the  only  requirement  to  he  fulfilled  in  the 
development  of  ocean  and  atmosphere  models.  At  the  present  stage  of  technology, 
computer  efficiency  might  he  an  even  more  restrictive  condition 

The  shallow  water  equations  are  the  prototype  equations  for  primitive  equation  nm 
dels  of  ocean  atmosphere  dynamics  It  is  well  known  that  numerical  techniques  ha 
se»l  on  explicit  time  differencing  schemes  aTe  considerably  affected  by  the  Cmirant 
Friedrichs  I  eyy  ((TL)  stability  condition  that  controls  the  high-frequency  gravity 
wave  motion  Although  explicit  schemes  are  always  more  accurate  than  implicit, 
t lie  latter  are  widely  applied  because  they  are  able  to  use  much  larger  time  steps 
(O’Brien,  I *)Xf» )  In  general.  fuUy-miphcit  schemes  are  seldom  applied,  and  semi 
implicit  schemes  are  often  applied  m  ocean  atmosphere  models  (Hamilton.  1977) 
Such  schemes  treat  the  terms  that  govern  the  fast  gravity  waves  miphcitJv  and  th* 
remainder  explicitly. 

The  implicit /semi  implicit  schemes  usually  require  the  solution  of  a  two  dimensional 
Poisson  or  Helmholtz  equation  A  large  number  of  different  solver  routines  exist  for 
such  an  elliptic  equation  However,  all  these  methods  are  time-consuming  and  most 
of  them  are  applicable  only  for  particular  boundary  shapes  or  boundary  condition 
specifications.  In  order  to  avoid  this  problem.  Tanguay  and  Robert  (19Kb)  have 
recently  proposed  an  algorithm  called  tin-  fractional  step  method,  which  reduces 
the  matrix  from  the  Helmholtz  equation  to  a  simple  tndiagmial  system  m  each 
of  the  two  horizontal  spatial  dimensions  Consequently,  solutions  can  he  efficiently 
obtained  by  a  special  form  of  ( iaiissian  elimination  met  hod  {( ’ar.iahan  et  al  .  I9t»q) 
The  practical  advantage  is  that  no  additional  computer  time  is  required  in  this  part 
of  the  calculations  (compared  to  the  total  computation  of  large  scale  dynamics 
models)  Although  the  languay- Robert  algorithm  is  based  on  a  modified  versi.in 
of  the  shallow  water  equations  with  the  inclusion  of  an  extra  nonphysical  term, 
the  method  looks  so  promising  that  the  main  purpose  of  this  study  is  to  verify  the 
accuracy  of  the  scheme  in  long-term  numerical  simulations  The  scheme  i«  therefore 
tested  and  compared  with  the  explicit  and  the  semi  implicit  methods 
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Section  2  presents  the  formulation  of  the  different  numerical  schemes.  The  schemes 
are  compared  in  Sect.  3,  and  the  results  discussed  in  Sect  4. 
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2.  The  scheme  formulations 


2.1.  THE  CONTINUOUS  AND  DISCRETE  SHALLOW  WATER  EQUATIONS 


We  consider  a  simplified  formulation  of  the  shallow  water  equations  (Pedlosky. 
1979).  The  equations  are  referred  on  a  /-plane  with  a  cartesian  coordinate  sy¬ 
stem  (x,  y)  chosen  such  that  in  the  Northern  hemisphere  the  ^-coordinate  increases 
eastwards  and  the  y-coordinate  increases  polewards,  taking  the  form 


1*1  /o«  -  SVx, 

(2.11a) 

>’(  +  /«“  =  -9'ho 

(2.1.lb) 

»,  d  t',)  -  0. 

(2  lie) 

The  subscripts  [?,y,1)  denote  partial  differentiation;  the  variables  (n,  v)  are  th»j 
components  of  the  eastward  and  poleward  velocities,  respectively.  The  variable  ?/ 
represents  the  free  surface  displacement,  g  is  the  gravitational  acceleration,  /o  the 
Coriolis  parameter,  and  H  the  total  depth  of  the  water  column.  Without  a  loss  of 
generality  we  assume  flat  bottom  topography  The  equations  are  satisfied  in  the 
domain  P  .  {(j\y)|  0  <  x  <  LT>  0  y  •  Ly] 


The  analytical  solution  of  Eq.  (2.1  1 )  may  have  the  simple  form 
(u.r.rj)  -  (f'.r,  N  "t*'"'. 


(2.1  J) 


where  JV  is  the  amplitude  of  the  wave,  uj  the  frequency,  and  k  and  (  the  wave 
numbers  It  can  be  easily  proved  that  the  coefficients  V  and  V'  as  a  function  of 
amplitude  are  given  by 


V  _ 


g  N 


v 


gN 

- 7j(^  *  "'/O) 


,2  '  mw  "jv/i  •  2  f2 

*  Jo  *  Jo 

The  wave  expressed  by  (2.1.2)  satisfies  the  dispersion  relationship 

-=  \Jfl  + 

and  has  a  group  velocity  GV  of  components 


GV 


_ _ gH±_  _  \ 

y/Ti y/ft+gHp)’ 


gff k 


(2.1.3a) 


(2.13b) 


(2.1.3c | 
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where  0  =  \/k2  f  /V  Only  positive  frequency  values  are  considered  henceforth 

The  analytical  problem  (2.1.1)  is  numerically  approximated  on  a  C'-grid.  We  have 
chosen  this  grid  because  it  is  the  one  most  often  used  in  large-scale  dynamical 
models  (Grammeltvedt,  1969).  All  the  terms  are  centred  differences  in  space,  with 
a  grid  network  such  that  the  functions  u,  t\  and  tj  have  the  numerical  correspondent 

u",m  =  U((‘~J  +  l)Ax,2mAy,nAf), 
t'"m  -f(2jAf,(2m  +  1  )Ay,  n Af ), 

'/j.m  —  »l(  2  j  Ar,2mAy,nAf ),  (2.1.4) 

where  j ,  »»t  and  ri  are  the  indices  relative  to  the  variables  t ,  y  and  f,  respectively. 
However,  we  will  suppress  the  indices  when  not  incremented.  With  the  above 
assumptions  and  notations,  the  solution  (2.1.1)  is  discretized  as  follows: 

U  =  Ar ),  ,,  _  yyt*.  rtayl,  _ 

<t>.  =  <t>"  „  =  ui.nAt  +  «2jAr  +  f2tnAy,  (2.1.5) 

where  *  indicates  values  to  be  computed  for  different  schemes.  The  schemes  are 
illustrated  in  the  following  sections 


2.2  THE  EXPLICIT  SCHEME 

With  the  explicit  scheme,  Eqs.  (2.1.1)  are  written  as  follows: 

untl  =  u"- 

1  +  2/oAflt’]"  -  ~(t;J  +  t  -  rJJ)n, 

(2.2.1a) 

fn+l  =  t-" 

1  -  2/oAf'  ul"  -  -(r/m  +  t  -  i/m  )n, 

Ay 

(2.2.1b) 

tr1  =  -r 

1  -  h Af( ,lj  Ui  '  + r"  1)~. 

Ax  Ay 

(22.1c) 

where  [.  ]  indicates  the  average  over  the  four  closest  points 

This  scheme  is  affected  by  the  CFL  stability  condition,  which  requires  a  tune  step 
Af  such  that 


V 
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Relative  In  lh*‘  solutions  (2  1.5),  the  scheme  gives  the  following  values: 


where 


«  V  (- 
El  A7  t  f 


n  sAr  ,  siii/Aji  \ 
tA 

Ac  Ay  ' 


r-  ("f 


Ay 


.  i  A 


Ax 


At 


'(Aty  A7  .  I'7) 


(2.2  3a) 

(2.2  31.1 


A  /<,  rns  a  Ax  ros  (Ay. 


1 


/  ,,  Sill 

V''‘  Ax 


A.r  sin*  /Ay 

*  Vs) 

Ay- 


sinwt-At 

At 


(2.2.3c) 


I  he  components  of  the  vector  group  velocity  are  symmetric,  so  we  will  consider 
only  one  of  them,  namely  the  eastward  ci  mponent  I.'nfort unately,  the  mathema¬ 
tical  representation  of  the  computational  group  velocities  might  he  quite  complex, 
and  although  intriguing  would  add  no  significant  information  about  the  associated 
numerical  distortions  Thus  it  is  only  the  group  velocity  values  for  non- rot  at  ional 
Hows  (viz  /n  H)  that  are  presented  henceforth 


<;vr. 


ijll  sm  2s  Ax 

l\(l  At7!'7)'  2  Ax 


•2.3d) 


We  recognize  a  formal  similarity  between  Kqs.  (17.2.3)  und  (2.1.3)  in  the  sense  there 
are  the  ‘computational  wavenumbers’  sin  kSt /  Ax  and  sin /Ay/ Ay  corresponding 
tn  the  ‘true  wavenumbers'  *  and  t ,  and  the  'computational  Coriolis  frequency*  A 
corresponding  to  the  ‘true  Coriolis  parameter’,  /o  These  computational  terms 
are  independent  of  the  explicit  treatment  of  the  Eq.  (2  11)  They  are  a  direct 
consequence  of  the  rentred  in  space  finite  difference  approximation,  and  the  (‘-grid 
(i  e  the  4'point  average),  respectively.  Besides  these  computational  variables,  there 
are  the  additional  computational  frequencies,  z.y.  and  a./.  While  the  frequency  Ey 
is  the  result  of  the  rent  red- in  space  finite  difference  scheme,  the  frequency  u,<E  is  in 
trmsir  the  explicit  formulation  Thus  we  takeu^/r  to  be  the  effective  ‘computational 
phase  frequency’  of  the  explicit  scheme 

It  is  easy  to  verify  that  all  the  computational  variables  converge  and  will  converge  to 
the  corresponding  analytical  values  as  the  increments  At.  A/,  and  Ay  tend  to  zero 


5 


SaC'LANTC'EN  SM-202 


It  generally  follows  that  the  representation  is  fairly  good  for  long  waves  (i.e  those 
resolved  by  several  grid  points),  but  poor  for  shorter  waves,  particularly  those  with 
wavelengths  less  than  4A.r,  4 Ay  a  wavelength  4  A*  corresponds  to  kAj*  -  jt/2 
(Grotjahn  and  O'Brien,  1976). 


J.3.  THE  SEMI-IMPLICIT  SCHEME 


The  semi-implicit  treatment  of  the  shallow  water  equations  is  written  as  follows: 


bfU" 

1 

K’l”)  -  /o1-''". 

dtl"  - 

fy’l’"' 

•  In'  i‘". 

- 

'  2  //(*,  u"~' 

■  dTu”  -  *-  <“!w i  ”). 

(2.3.1) 

where  <*,,  and  bs.  are  the  centered  finite  difference  operators  in  time  and  space 


respectively.  With  the  use  of  incremented  index  variables,  Eq. 

(23.1) 

are  written 

as 

-  1j)n" 

=  01. 

(2.3.2a) 

02, 

(2.3.2b) 

n  ,  11  -if 

1  //Af 

--  03. 

'1  -1 

)  •  -  .  ( i1...  -  - ,  ) 

2Ay 

(2.3.2c) 

where  Q 1,  Q  2,  and  Q  3  contain  the  remaining  terms  evaluated  at  time  step  nor  ri-  1. 
Substitution  of  (2.3  2a)-(2.3.2b)  into  (2.3.2c)  leads  to  the  numerical  Helmholtz 
equation  for  7/n*1,  when  all  the  variables  at  time  step  n  are  known: 


-a’lj+l  O'h-l  -  H n.i  b’lm  1  -»  (l  *  -<>  :  -M>;  -  Q. 


(2. 3. 2d) 


(2.3.2e) 


where 

y//A/2  ,  qH&t2 

a  ‘ — — -  ,  b  ■  ' - •-  . 

4A-r*  4Ay* 

Once  the  surface  elevations  are  determined,  the  velocity  field  can  be  updated  from 
Eqs.  (2.3.2a)-{ 2.3.2b).  The  scheme  is  still  affected  by  the  CFL  stability  condition 


/0JW  •  1 


With  respect  to  the  analytical  solution  of  (2.1  1).  the  scheme  gives  the  following 
computational  values. 


y:V  cos(i*;S/Af)  /  _  sin  k&i 
H2/  A2  \  Sl  Aj* 


I A 


sin  (  Ay 


S.A<  I  ANT<  'FN  SM.:»02 


./A  '  ;  sm'Aw  sin  kA/ 

..  z.  <  /  •  /  A  A 

E*,  A'  1  Ay  A.r 


I  A?  •  I  - 

■  s/  „  «trt-»  an  A/ «  4  „ 

A/  \  1  A/:.\- 


rj  :t  :u. 


'A,/ 


'/  H  vin  L'kA/ 
•  A/:Pl  2Ar 


CJ  :t  :t* 


uli-’r*  .r - 


;s  I  Tinalh  id*t:»  j'  al  t •  >  z  / 


A  formal  comparison  lw>t\vn*»n  I>jv  (23  A)  and  (2.1  il)  indicates  (hat  tlm  ■«*,nii 
implicit  sclmm*4  affects  t li«*  amplitude  of  tin-  wav**,  introducing  also  a  'computa¬ 
tional  amplit  u « 1  *r*  ’  A  .•  A' ’  <  oM  ,*«,/ AM  Thus.  '  lm  s^nn  implicit  sclmum  t*»nds  to 
mi  *  l*'r*^st  tmat*-  i>«th  l"r*"ip  i*  ’ii  c  \  and  ampin  .  i  *  I  *  ■  of  tit*-  waves 


;  !  It!  I  HA*  I  i.  >'•  M.  'll  i*  Ml  1  IP  ‘I  ■ 

In  -.pf*T  la-  ton;*  tin-  H-  imliolt/  1  |  i  2  2a  t.  tin  fra*  ti«>nal  slop  method  suits*  i 
t  nt*  -  f  It-  •  <  ■miiiun  v  I  j  *2  i  !  •  I  wit  li 

•  n  ■  \J  i .  •  II  [  n ,  •  v,t  i  0.  (2.11a) 

vvlnr* 

-  •»  *<  *  u 

M  I  .///A/*  r  „  ,  (2.4  lh) 

lr*atimi  tin-  lap  '2  I  la.h*.  ami  <2.4  li  '“tni-implicit  lv.  I  In4  rporpiatinn  (2.12a) 
t  t  In-  !•  Tin 

L  i  r;  i  L'yL'^h])  (1  •jHAt't'l  |(1  </  //  A/2**  J )tj  Q  + pM'1.  (2.4  2) 

w  in  r*-  C  r .  and  C u  ar*  tin-  «;*•'•« »n< I  ord**r  **-nl«-r*,d  difference  operators  mi  -  and  y. 
r<  sp-  .-t  i v # •  I v .  ami  p  ( */ //  A/:  7/  Since  tin-  t«*rni  p  is  « »f  fourth  ord*»r  in  A/, 

tin  srln  mi"  sln-uld  not  introduce  aiiv  Muniln  ant  '-rr*»r  in  tin-  solution  (  fan^uav 
an*l  K*.l**Tt .  I HKti  i 


:>ut  at  tonal  vain*'-  associated  to  t  h*  srlu'iiu-  an 


<J  V  r,is(^f  A/ ) 

=\  \- 

y  A  cosf  u.’/.  A  / ) 


'in  k  A  j- 


»A 


sin  ( A  t/ 


rnsia,.  at)  / 

=i-  a*  ' " 


A  j*  A  y  ' 

sin  (  At/  sin  k  A  j  ' 


Ay 


i  A 


A.r 


'2  I  21a  | 
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r 


At 


arctan  At 


A  =  (  I  •  At2'i2 )  •  I  - 
\  (I  •  At-n'-H  !  At'-  A:)  ! 


( 21  :«> 


<  ;v 


1  •  At-'i-  i 
r  1  .  .  P|  \  f  2 A r 


y  H  sin  2  k  A  x 


I'o  sin  /Ay  z 


wher 


yffAt 


siiikAj*  Mil  /  Ai/ 


A-r 


Ay 


Af\y  //  <■< iskA/  j  . 

(2  l  :ir| 

(2.1  ;Mi 


It  is  easy  to  verify  that  the  fractional  step  computational  variables  converge  to 
the  semi-implicit  computational  variables  as  ()( At4  )  Therefore  for  small  At  incre¬ 
ments  the  two  schemes  should  end  up  virtually  identical,  ronfirming  the  Tanguav 
.and  Robert  (1980)  hypothesis. 

The  r-  ami  y-sviniuet  rv  of  the  operator  C  of  (2.4  2)  suggests  the  solution  of  the 
rj-equation  using  an  alternating  direction  iterative  algorithm  in  which  the  x-  and 
y-operators  are  inverted  at  each  time  <uep.  At  each  iteration  only  the  boundary 
conditions  at  two  opposite  sides  are  necessary.  Let  us  solve  the  operator  (2.4.2) 
with  the  following  algorithm 

First,  we  compute  from 

CTv  (I  ytfA/‘^)v  Q  *  :  (2.4.4a) 

at  the  inner  ?/  mesh  points;  no  boundary  conditions  are  required  at  this  stage  except 
for  the  terms  that  appear  in  the  righhand  side  of  (2.4.2).  Now  we  compute 
from 

(1  gHM2il)n  v  (2.4.4b) 

using  only  the  boundary  conditions  at  the  northern  and  southern  boundaries.  At 
the  following  time  step  the  operator  CT,  and  Cy  are  inverted.  Subsequently,  the 
use  of  an  alternating  direction  algorithm  allows  a  more  accurate  application  of  all 
the  boundary  conditions. 


R 
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3.  Applications 


To  verify  the  accuracy  of  the  schemes  fully  we  now  compare  numerical  results  with 
exact  analytical  solutions  of  problem  (2.1.1 ).  Numerical  experiments  are  performed 
with  the  following  common  features:  total  depth  H  —  2000  in,  total  extensions 
/. j  1000  km,  and  Lv  ■  1000  km,  and  grill  spacing  Ac  17  4  km,  Ay  . 
17  4  km  Although  the  formulas  presented  in  the  previous  sections  are  formally 
symmetric  with  respect  to  the  wavenumbers  *  and  t,  there  is  a  different  structure 
to  the  errors  because  they  contain  quadratic  terms.  Nevertheless,  the  computational 
frequencies  and  group  velocities  of  the  schemes  are  depicted  in  a  one-dimensional 
form  in  Figs.  I  3.  For  simplicity  we  consider  irrotational  flows  (ie.  /o  0)  and 
choose  the  wavenumbers  such  that  sAr  fAy  The  schemes  are  tested  for  various 
values  of  the  non-dimensionalized  ratio 


-  VyH 


At 

v'2Ar 


Cirotjahn  and  O'Brien  (1976)  gave  an  accurate  and  complete  analysis  of  the  di¬ 
stortions  introduced  by  the  explicit  and  semi-implicit  time  differencing  schemes. 
In  general,  explicit  formulations  have  the  tendency  to  overestimate  the  oscilla¬ 
tions.  whereas  the  implicit  have  the  tendency  to  underestimate.  Although  explicit 
schemes  are  found  to  be  more  accurate  than  the  implicit,  the  accuracy  is  usually 
incommensurate  with  the  higher  computational  cost.  The  same  findings  apply  to 
the  fractional  step  method  As  Fig  3  confirms,  the  fractional  step  method  intro¬ 
duces  the  same  distortions  of  the  semi-implicit  approximation,  but  it  is  slightly 
less  accurate.  However  the  errors  might  well  be  compensated  for  the  increased 
computational  efficiency. 

Finally  we  compare  the  solutions  obtained  bv  solving  numerically  the  problem 
(2.1.1 )  with  the  exact  solution 

(«.  <’,!))  (F,  I’,  A  )sin(wof  t  Kol-r  >  y)), 

where  the  wavenumber  is  k0  7  86  x  If)  8  cm  1 ,  and  the  frequency  is  w0 
1 .44  s  10"*  s'  1 .  The  amplitude  is  iV  -  20  cm,  the  velocities  V  and  V  are  defined  as 
in  (2.1  3a)  The  wavenumber  k0  corresponds  to  a  wavelength  A  ---  50  in  terms  of  the 
grid  intervals  To  obtain  a  good  resolution  of  the  wave  period,  the  implicit  schemes 
used  a  time  step  At  — -  217  s,  which  implies  20  time-iterations' per  period,  and  a 
value  T  —  1 .25  The  explicit  scheme  employ  a  time  step  At  --  50  s,  corresponding  to 
a  value  r  ~-0.3. 
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Fig  1:  The  computational  frequency  (a)  and  com}  (na¬ 
tional  group  velocity  (b)  of  the  explicit  scheme  norma¬ 
lized  by  the  analytical  exact  value.  The  abscissa  repre 
sents  the  variable  kAx  normalized  by  jr.  (I)  r  '1.0 1, 
or  At  -  1.75;  (II)  r  -  0.5  or  Af  --  87. 
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DO  (il  C  5  0  6  1.0 


/■*! jl.*  J  lli*-  >  >m (>ui  i< >ii.il  I r | »i *- ii *  \  I  <t  \  4 ii* i  < 
h-'ii.tl  t  r  •  •  ii  |  •  Hi  i  hi  <1  t  Ii**  »»nu  implicit  s«  li**tn*  jl> 

r  M  *>  t  .if  A/  IT',  till  r  II  ..r  HT.  (IMl 

-  I  •*[  ±1  I  ”  •  i  I  \  ■  *■  in  **r  I  "Ml 


Group  Velocity 
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0  0  0.3  0.6  O.e  1.0 

F’jg  J  The  computational  frequency  (a),  and  computa¬ 
tional  group  velocity  (b)  of  the  fractional  step  method 
Same  a s  Fig  l. 
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Ki^ur**  illustrates  fh*  surface  •  1 1 *» t rdnitp >n  ■■\«r  tie  whole  basin  :il t <*r 

*  wave  p*»n<"N  It  i"  evident  that  **{»♦*<•  i ft*- ?it  i* »n  of  ih»*  !•«  h i ii< ! -ir v  •' '*  »ndit  i"ii>  together 
with  the  distortion  between  computational  frequency  1 1  ♦*  t  Ip- frequency  “I  tip- mi* 
nor)  and  exact  frequency  ( i  ••  tin*  frequency  at  tip-  boundaries)  afb'ts  tin-  quality 

•  »f  t |i«*  solution  In  ip>rai  all  the  numerical  schemes  ar*  sensible  to  this  prol.|.  in 
High  frequency-  waves  that  aro  present  at  tit*-  boundaries.  mu;iit  l>.  remove.!  fa 

I  transforming  the  boundary  conditions  into  a  lore  mi;  term  for  tip  r/  .-quaiion  I  his 

could  he  achieve*!  by  introducing  a  new  function  ?/  where  */  is  an 

arbitrarily  chosen  high  frequency  function  such  that  »/'  H  on  .</»,  and  solving  the 
differential  problem  for  if'  In  this  a  way  tlv  solution  is  locally  adjusted  to  tip- 
short  waves,  and  the  (listort ions  between  computational  and  exact  frequences  ar* 
notably  reduced 
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